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The nonlinear dynamics of preheating after early-Universe inflation is often studied with lattice 
simulations. In this work I present a new lattice code HLATTICE. It differs from previous public 
available codes in the following three aspects: (i) A much higher accuracy is achieved with a modified 
sixth-order symplectic integrator; (ii) scalar, vector, and tensor metric perturbations in synchronous 
gauge and their feedback to the dynamics of scalar fields are all included; (iii) the code uses a 
projector that completely removes the scalar and vector components defined by the discrete spatial 
derivatives. Such a generic code can have wide range of applications. As an example, gravity waves 
from preheating after inflation are calculated with a better accuracy. 



I. INTRODUCTION 

Early-Universe inflation |l|-l3| has become one of the 
key elements of the standard cosmological model 
In this paradigm the Universe went through a phase of 
accelerated expansion, which, in the simplest scenario is, 
driven by a scalar field, namely the inflaton. The pre- 
dictions of inflation have been confirmed by the high- 
precision measurements of CMB 0413] and the large 
scale structure surveys jlS^ffl. Combining these obser- 
vations with supernova [211428]. wea k gr avitational lens- 
ing (29rl36| and Lyman-a forest data [37H4l] j , we find that 
the current Universe is undergoing another cosmic ac- 
celeratio n |42T - (44| . This could be due to a cosmological 
constant \4M or another scalar field |46to7|. 



It is hence important to understand the dynamics 
of scalar field(s) in a perturbed Friedmann-Robertson- 
Walker (FRW) background In many cases these 

scalar fields are almost homogeneous. Thus linear or 
second-order perturbation theory is enough to describe 
such a system. However, there are exceptions. Some- 
thing that we cannot avoid in any successful inflationary 
model is the decay of the inflaton after inflation. This 
could start with preheating, a nonperturbative violent 
process due to parametric resonance or tachyonic growth 
of fluctuations [58h69| • The typical comoving scale of pre- 
heating is much smaller than current cosmological scales. 
However, cosmological-scale comoving curvature fluctua- 
tions can also be generated via, e.g., preheating mod- 
ulated by a field that is light durin g in flation but be- 
comes heavy at the end of inflation [70l - l72l |. Moreover, 
the stochastic background of gravity waves (GWs) gen- 
erated during preheating [73148 ll | is, in principle, observ- 
able. In particular, GWs from preheating after hybrid 
inflation may be observable with the next generation of 
GW probes (79|, although this depends on the parame- 
ters in this model. See also (75[, where the parameter 
space is systematically explored. For more models sug- 
gesting potential ob servab les from preheating, the reader 
is referred to Refs. 



To make quantitative predictions of the observables 
from preheating, one often needs to run full nonlinear 



lattice simulations. In the previously mentioned GW cal- 
culations, the evolution of scalar fields is done in con- 
figuration space using the public available code LAT- 
TICEEASY [88[ or unreleased codes with similar tech- 
niques. The linear metric perturbations are evolved ei- 
ther in configuration space (75l - [79| or in Fourier space 
[z3, HSIHI- In these treatments, a traceless-transverse 
(TT) "tensor mode" is defined in Fourier space 



iTT 

n ij,k 



[M TT (k)l... h, 



(i) 



where hij are the metric perturbations in synchronous 
gauge 0, |90[ • The Fourier-space matrix form of the TT 
projector, M TT (k), is given by 



[M TT (k)] .. lm ee P a (k) P mj (k) - (k) P lm (k) 



where 



Pi 



00 - % - ^ 



(2) 



(3) 



Here 8ij is the Kronecker delta. If not otherwise stated, 
repeated indices are implicitly summed over. The Latin 
indices run from 1 to 3 (spatial indices). The Greek in- 
dices run from to 3 (temporal and spatial indices). 

One could also define the TT component in configura- 
tion space: 
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where 



V 2 = dl + d\ 



h = hi. 



diAi = 0, dihtf = 



(4) 

(5) 
(0) 



Here hij is decomposed into two scalar components (h 
and A), one vector component {Ai) and one tensor com- 
ponent (hJj T ). The tracelessness of hj^ can be confirmed 
by taking the trace of Eq. (0| and by using Eq. ^ . 

In the continuous case, definition ((?]) is equivalent to 
definition |T]). But this is not so for a periodic and cubical 
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box with length L and n 3 grid points, in which the spa- 
tial derivatives in Eq. (|]D) are replaced by finite difference. 
This is discussed in detail in Sec. Ill Bl The discrepancy 
between definitions ([T]) and @ leads to scalar-tensor mix- 
ing: the scalar part of the energy-momentum tensor, cal- 
culated in configuration space using discrete derivatives 
(finite difference), can produce GW in Fourier space de- 
fined by Eq. ([1]). At scales where the scalar components 
dominate, it is difficult to suppress this "noise" or dis- 
tinguish it from the physical GW. Simulations with very 
high resolution (n > fO 4 ), which might solve the prob- 
lem, are not favored, practically, as they are numerically 
expensive and often limited by the machine memory. In 
the new code HLATTICE that I will present in the paper, 
I define, evolve, and extract scalar, vector, and tensor 
modes consistently in configuration space. Even though 
the discrepancy between real (continuous) physics and 
the numerical (discrete) model cannot be completely re- 
moved in any numerical calculations, in HLATTICE the 
scalar, vector and tensor parts of the metric perturba- 
tions are only sourced, respectively, by the scalar, vector, 
and tensor parts of the energy-momentum tensor. 

In other public available lattice codes, such as LAT- 
TICEEASY, DEFROST [U, and CUDAEasy the 
scalar fields are evolved in a FRW background, and met- 
ric perturbations are ignored. But at linear level one can 
approximately take the energy-momentum tensor of the 
scalar fields as a source, and evolve GW outside the sim- 
ulation. (Therefore, the scalar-tensor mixing effect is a 
problem in post-processing. It should not be regarded as 
a problem of these lattice codes.) Since this is a linear 
treatment, the TT component separation could be done 
at the end of the calculation [79(. Except for Ref. [93| . 
which I will discuss separately in Sec. IIV1 the previous 
works on GWs from preheating are all based on this (or 
a similar) approach. In these calculations, the feedback 
of metric perturbations to the dynamics of scalar fields is 
ignored or partially ignored. HLATTICE is the first code 
released that consistently evolve all components of met- 
ric perturbations together with the scalar fields. Using 
HLATTICE I find the metric feedback, as conjectured 
in previous works, is indeed not a dominating effect, at 
least for the models studied in this paper. 

In some situations, in order to capture some small ef- 
fects [72| or energy- insensitive modes [94 |. we need to 
evolve the system of scalar fields and metric perturba- 
tions accurately. In HLATTICE I use a sixth-order sym- 
plectic integrator for global evolution, and a fourth-order 
Runge-Kutta integrator 95] with refined time steps for 
the noncanonical terms only. The advantage of doing 
so is that no extra memory cost is required. With this 
integrator the fractional energy noise of the system can 
be suppressed to < 10 -12 . This enables us to check the 
conservation of the total Hamiltonian, including the tiny 
contribution from energy carried by gravity waves. This 
is the first time that we can use the constraint equation to 
accurately check the numerical accuracy in calculations 
of GW from preheating. 



This paper is organized as follows. In Sec. |TT] I intro- 
duce the HLATTICE code; in Sec. HID I use HLATTICE 
to calculate GW from preheating. I discuss and conclude 
in Sec. HVl 



II. HLATTICE CODE 
A. Theory 

The system that we consider contains m canon- 
ical scalar fields <j>\, (j>2, (f> m with a potential 
V (<t>\, <t>2, 4>m)- The action reads 



S = J ^d 4 x(^ 



Ml 



d^d^t-V+^R] , (7) 



where g = |det <? M „|, g^ v being the spacetime metric; M p 
is the reduced Planck mass, related to Newton's gravita- 
tional constant Gn via M p = 1/\/8itGn', R is the Ricci 
scalar Q. 

The _sp_acetime metric can be written in synchronous 
gauge 



ds 2 = dt 2 — gijdx l dx^ 



(8) 



Natural units c — h = 1 are used. 

In the context of inflation, we are interested in an ex- 
panding and spatially flat Universe, where the metric gij 
is often written as 



gij = a(t) 2 (Sij + hij) 



(9) 



In a finite volume L 3 , we choose a(t) to be the "scale 
factor," given by 



<*(*) = (Z3 / 



1/3 



The Hubble parameter is defined as 



H =- , 
a 



(10) 



(11) 



where a dot denotes the derivative with respect to time 
coordinate t. 

We will consider weak gravitational fields that satisfy 
the condition \hij\ <C 1. For a problem that requires 
long-time integration, a growing nonphysical gauge mode 
might spoil this condition. This specific case is discussed 
in Sec. lUDl 

In HLATTICE I choose to evolve the following vari- 
ables 



(12) 



where Q is the 3x3 matrix g^ . The matrix function In Q 
should be understood as 

InS = 2£\na+{a- 2 g-£)-\{a- 2 g-E) 2 +\{a- 2 g-Ef-. 

(13) 
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where £ is the 3x3 identity matrix. To the linear order, 
we have 



fcj - 25ij In a 



(14) 



It follows from Eq. ([14)) that the traceless part of fiij, 
which we define as 7^, satisfies 



V 



h 

3' 



Hj = Pij ~ = hij - -5 l3 + 0{h%) . (15) 



This implies that 0(7™) < 0(h™) for arbitrary n > 0. 
This significantly simplifies the procedure to expand an 
arbitrary function of fiij to a given order in hij- we just 
need to replace f3ij with ^Sij + jij and cut the Taylor 
expansions in 7^ at the same order. In the rest of this 
subsection, we will do this for the metric and for the 
action. 

Let us first calculate the volume weight Jg. Using the 
fact that dct(Cf) = e Tv ^ nQ \ we find a simple and exact 
expression: 



V9 



- p/3/2 



(16) 



where j3 = 0a is the trace of /3y . 

Writing the 3x3 matrix 7^ as T, we can rewrite the 
3x3 matrix <?jj as 



g = Je+T = e f>/3 ^ + p + lp2 + lp3 

By taking the inverse of Q we obtain 



(17) 



9 ij = e~° Uii ~ Hi + f <<■"<.•<) + O (III) ■ (18) 

We can substitute, for example, 1 — 711 + ^Tii with 
e -7n _|_ O(h^). With such substitutions and ignoring 
0(ft| ) terms, we can write g y as functions of 



,11 w e -ft 



22 w e -/9 22 + 



,,33 



e -2^u/3 
2 

e -2/3 22 /3 

2 

e -2& 3 /3 



}2 p -/3 22 /3 



7 23 « -/3 23 e-^-+fe)/ 2 + i^^aie-^ 3 
? 12 « -ft 2 e-(/ 3 -+fe)/2 + i^ifee^ 73 



(19) 



The difference between the left-hand side and the right- 
hand side of each equation is less than or equal to 0(h^). 
It is obvious that such approximations are not unique, as 
one can add arbitrarily more higher-order terms on the 
right-hand side of each equation. The specific choice I 
made in Eqs. (|19p is based on three criteria: 



• simplicity; 

• less 0(/i 3 ,) residual error terms; 



• under a coordinate transformation x % — ¥ Cx l (i = 
1, 2, 3, C is a constant; the metric transforms as 
fiij — ¥ j3ij — 2Sij InC), the approximated expression 
of 3 U scales as C 2 , just as the exact <f J should. 

The last requirement enforces the gradient energy density 
g v di<f>gdj<f>g, with the approximations of g^ in Eqs. (fl"9|) 
being used, to remain exactly invariant. This allows us 
to perform a spatial coordinate transformation x % — > Cx l 
without producing any extra error terms. This can be 
used to optimize the HLATTICE code. In HLATTICE, 
after every few evolution steps the coordinate redefini- 
tion x l — > Cx % is performed, where C is chosen to be 
the scale factor at the moment. After such a transforma- 
tion, the scale factor is redefined to be 1 and fiij is made 
small. The exponential functions in Eqs. (I19|) can thus 
be evaluated via the approximation 




(i/16) s 



for \x\ < 1. (20) 



For a program built with Fortran code, the evaluation 
of the right-hand side of (f2"0)) is about twice as fast as 
that of e x . This significantly improves the performance 
of HLATTICE. 

Confusion should be avoided here. The FRW back- 
ground is intrinsically different from a Minkowski back- 
ground. In the calculation, by adaptively changing the 
coordinate system, it is possible to keep the scale factor 
close to 1 and to keep /3ij small. However, when deriv- 
ing the theoretical formulas, we should work in a fixed 
coordinate system and cannot assume that /Jy is small. 

The action ([7]) can be written as 

S = Jdt(K f -G f -V f + K g - G g ) , (21) 

where Kf is the kinetic energy of the scalar fields, 

1 



K f = / e p/z d J x 



2j3„ ± J2 



Gf is the gradient energy of the scalar fields, 



(22) 



G f = J eP/Wx^difadjfa ; (23) 
Vf is the potential energy of the scalar fields, 

V f = J Wi^,^ W ; (24) 

Kg is the "kinetic energy" of gravity, approximated to 
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second order in hij by 
M 2 r 

X ($3 + 031 +012 

-011022 - $22033 ~ 0330n) ! (25) 

and G g is the "gradient energy" of gravity, approximated 
to second order in hij by 

G g pb -^-o(*) / d 3 z 

X (023,1 + 031,2 + 012,3 

— 2^23,1(331,2 — 2p3i,2012,3 — 2/?i2, 3023,1 
~022, 1033,1 — 033,2011,2 — 011,3022,3 
+2/323,2011,3 + 2/331,3/322,1 + 2/?i2, 1033,2) , (26) 

where pV/,fc = dk(3ij- In Eq. ([26|) I have approximated the 
local volume weight e^/ 6 with a global scale factor a(t). 
This is valid since the integrand in Eq. (|26p is of sec- 
ond order in hij. In Eq. (|25l) such a replacement is not 
allowed, as the integrand is of zeroth order. With the 
bilinear approximation (|2^|) . HLATTICE cannot capture 
the gravity self-interaction, which in principle can be im- 
portant on small scales. This is discussed in Sec. IIV1 

Equations (|25M26p are obtained using the same tech- 
nique that I used to derive Eqs. (fll)]) . i.e., rewriting /3y 
as jij + ^5ij and cutting the Taylor expansions of 7^ 
at the second order. The simplicity of the final expres- 
sion of Kg is the main reason why I have used pV, as 
fundamental variables in HLATTICE. 

It is also easy to verify that the right-hand sides of 
Eqs. (|22H26I) are all exactly invariant under the spatial- 
coordinate redefinition x % — > Cx % . 



B. The discretization scheme 



where n is a fixed integer representing the resolution of 
the simulation; / represents all physical quantities in- 
cluding the scalar fields, the metric, and their tempo- 
ral/spatial derivatives. Given the PBC, we only need to 
evolve fi 1 ,i 2 ,i 3 in a cubical fundamental box containing 
n 3 grid points (— n/2 < 11,12,13 < rl /2)- In what follows 
I will use the notation ^lattice ^° denote the summation 
over grid points in this fundamental box. 
The DFT is defined as [H 

7. . . = \ ~* p-^-iiiji+izh+iste) f . . . roe\ 

lattice 

where i is the imaginary unit and j± , j% , j'3 are arbitrary 
integers. Unless otherwise stated, the overhead tilde sign 
represents variables in Fourier space. 

It is easy to verify that in Fourier space fj lt j 2 ,j 3 also 
satisfies the PBC 

fjl+n,j 2 ,33 = fji,h+n,j 3 — /jij2j3+n = fji,h,h ■ (29) 

Moreover, if a field satisfies the PBC in Fourier space, 
it also satisfies the PBC in configuration space. Thus in 
what follows we no longer distinguish between the two 
PBCs. Similarly we define a Fourier-space fundamental 
box with -n/2 < ji,fa,ja < n/2. 

The standard DFT wave vector is defined as 

2"7T 

K:,,,. T Uuh,h), (30) 

where L is the coordinate length of each edge of the fun- 
damental box in configuration space. The amplitude of 

k-std •„ 

I will use the subscripts "disc" and "cont" to distin- 
guish between quantities evaluated on the discrete lattice 
or in a continuum. These subscripts are often omitted 
when it docs not give rise to ambiguity. 



In this subsection I answer or attempt to answer the 
following questions: 

1. What exactly is calculated in a lattice simulation? 

2. How do we choose a discretization scheme to in- 
clude metric perturbations? 

3. How do we define the scalar, vector, and tensor on 
the lattice? 

Before answering these questions, let us briefly review 
the lattice theory and the discrete Fourier transformation 
(DFT). 

I label a grid point in the lattice with three integer 
numbers (ii, i<i, %■&) and apply the periodic boundary con- 
dition (PBC) 



fii 



+n,i2, *3 



fii 



il+n,i 3 



fii 



i?M+n 



(27) 



1. What exactly is calculated in a lattice simulation? 

In LATTICEEASY, DEFROST, and CUDAEasy, the 
metric perturbations are ignored. Scalar fields are 
evolved in configuration space. The lattice- version equa- 
tion of motion (EOM) of the scalar fields is 



dt 2 



3H 



dV 



0, (32) 



where V 2 is the discrete Laplacian operator. 

In this subsection we will discuss an idealized case 
where the temporal differential operator d/dt can be per- 
fectly integrated. Moreover, the numerical errors in a(t) 
and other background quantities are assumed to be negli- 
gible. These assumptions can be quite close to the reality, 
if a high-order integrator is used in the lattice code. 
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The only spatial operator that appears in the EOM is 
the Laplacian operator V 2 , which can be directly defined 
without referring to any first-order discrete derivatives. 
For example, in LATTICEEASY V 2 is defined as 



V 2 / 



^2 (fii+hhM + fh-ijiM (33) 

fh,ji — l,fei 
j'i,fei + l "F fii,ji,ki—l ~ Qfii,ji,ki) ? 



where A = L/n is the coordinate distance between neigh- 
boring points in configuration space. By taking the DFT 
of the above equation, we obtain its equivalent form in 
Fourier space, 



V 2 / 



- - (k oS 
\ 31,32,33 



) In 



32,33 1 



(34) 



where 



/'/off \ z _ ■* / • 2 ^Ji i • 2 "^2 . 2 "J3 
jWa) = -X2 [ 8111 — + 8111 — + 8111 — 



2 TJ2 



,2 ^3 

n 



Because both V 2 / and / satisfy the PBC, their ratio 
-(fc cff ) 2 should also. This can be explicitly checked in 
Eq. (011). 

The reader should keep in mind that Eq. (j35j) is a 
consequence of the LATTICEEASY definition of discrete 
V 2 . In other lattice codes the definition of discrete V 2 
(and hence k eS ) can be different. 

Using Eq. we rewrite the EOM, Eq. ([32) m Fourier 
space: 



dt 2 



ueS 

31 ,32 ,33 



3H 



dV 



dt 



31,32,33 



. (36) 



31,32,33 



Equation (j3"6"|) is the equation being integrated on the 
lattice. Comparing it to the continuous-case EOM for 
a mode k = kjf ,- 2 ,- 3 , we find that the only discrepancy 
comes from the difference between the DFT of dV/dcj)e 
and the continuous Fourier transformation of it. This 
discrepancy is model dependent and difficult to quantify 
in the nonlinear regime. However, in the linear regime, 
as I will show below, this discrepancy vanishes. 

In the linear regime, we can rewrite the lattice- version 
EOM in configuration space as 

d 2 V 2 d 



dt 2 



a 

d 2 v 



dt 



■> *M llia ,i, = . (37) 

where (•) represents the lattice average ^ lattice '• The 
field perturbation 5(j>i is defined as 



(38) 



The Fourier-space version of Eq. (|37j) is 



dt 2 



31 ,32 ,33 



dt 



d 2 V 



31,32,33 



, (39) 



which is identical to the continuous-case EOM for a linear 
mode with wavenumber k — k cS . If we interpret k eS 
as the physical wavenumber, in the limit that a perfect 
integrator is used, all the linear modes will be correctly 
solved on the lattice. In other words, in a lattice code 
with a good integrator, the difference between a discrete 
system and a continuum does not matter in the linear 
regime. 

What has been discussed above is well known to the 
lattice community. However, in many previous works, 
£std m f unc l ain ental box is interpreted as the phys- 
ical wavenumber [74T - [8ll HI, IHj]. This is usually not a 
serious problem. With a good definition of discrete V 2 , 
k eS should be close to A: std in the fundamental box. How- 
ever, there are exceptions where the difference between 
k eS and A: std is relevant. One specific case is the calcula- 
tion of GWs that I will discuss later. 



2. How do we choose a discretization scheme to include 
metric perturbations? 

We have seen that defining the discrete V 2 in config- 
uration space is equivalent to defining a kernel — (fc ) 2 
in Fourier space. First of all, the kernel should satisfy 
the PBC. Second, since the linear modes on the lattice 
behave like linear modes with k = k eS in a continuum, 
we should restrict (k eS ) 2 to be positive to avoid non- 
physical tachyonic growth. Third, for practical purpose, 
in configuration space the discrete operator V 2 should 
involve only a few neighbors. Thus the functional form 
of (fc eff A) 2 is limited to be a low-order polynomial of 
e ±2i 7 rii/n ) e ±2i7rj 2 /n and e ±2i7rj 3 /n Finally, we want to 

choose a (k cS ) 2 that is close to (A: std ) 2 , as we want the 
nonlinear-regime mode-mode coupling on the lattice to 
mimic the real physics in a continuum. 

Using Eq. (f3"3"|) , the reader can verify that the LAT- 
TICEEASY discrete V 2 satisfies all the above conditions. 
Moreover, in configuration space it is accurate to the lin- 
ear order of A: 



V 2 /L = V 2 / 

d I disc d l 



0(A 2 ) 



(40) 



While this discretization scheme is simple and reason- 
ably accurate, there may be various reasons to improve it. 
For example, the discretization scheme in DEFROST im- 
proves the isotropy of V 2 without much additional com- 
putational cost [91|. In DEFROST the default discrete 
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V 2 is defined as 



where the effective wave vector ^jfj 2 j 3 is 



V fi lt i 2 ,i 3 — y ] C|j^ | + |»£| + |ij| fi 1 +i' 1 ,i 2 +i' 2 ,i3+i' 3 ' 

(41) 

where C = -64/15, Ci = 7/15, C* 2 = 1/10, and C 3 = 
1/30. 

For HLATTICE, which includes metric perturbations, 
however, the discretization problem is much more com- 
plicated, for the reasons that I list below. 

1. The governing equations for the metric perturba- 
tions are rather complicated. It is possible that 
nonphysical tachyonic growth arises in metric per- 
turbations even when the discrete V 2 is negative 
definite. 

2. The field EOMs now involve first-order derivatives 
of the fields, whose discrete forms should be prop- 
erly defined and should be consistent with the lat- 
tice V 2 . 

3. With the metric perturbations, the gradient en- 
ergy terms (G/ and G g ) are much more compli- 
cated. This prevent us from choosing complicated 
discretization schemes, whose computational cost 
would be intolerable. 

To address the first point, I run HLATTICE for empty 
spacetime with initial small noises in the metric, and ver- 
ify that the noises do not grow. I propose that such a 
null test should be done for any discretization schemes of 
gravity. 

Following the last two points, in the first released ver- 
sion of HLATTICE (HLATTICE VI. 0) I have used the 
simplest discrete V that is accurate to the linear order of 
A: 

^lli 1 ,i 2 ,i3 / = 2A ^' 1 +' 1 '* 2 '* 3 _ fh — ltiztis) > 



^ 2 /lil, 12,^3 



1 

2A 

1 

2A 



1,*2+1,*3 /«1,«2— 1,*3 ) , (42) 



{fn, 



fix, 



The equivalent Fourier-space form can be achieved by 
taking DFT of Eqs. @2J|. The result is: 



3i,32,]a \ \ .. / J 3i ,32, 3a > 
I . ( 27TJ2 



Kf ■ ■ ■ 

iJ 3i, 32,3a 
WJ J 31,32,33 



. ■ sin . 
A V n 



(43) 



Sm I I f 31, 32,33 



A 



The above equations can be written in a more compact 
form 



ik 



eff 

■31,32,33-131,32,33 > 



fix,. 



(44) 



,cir 



1 / . 2nj 1 . 2nj 2 . 2nj 3 

— sin , sin , sin 

A V n n n 



(45) 



Thus the Fourier-space kernel for —V 2 is 



i > 1 ( ■ 2 2-,/ 1 

) = ^2 ( S1E1 — 



27TJ2 



27Tj 3 



(46) 

Either taking the inverse DFT of Eq. (|4"6"|) or repeatedly 
using Eq. (j42p . we obtain the configuration-space defini- 
tion of discrete V 2 in HLATTICE VI. 0: 



V 2 / 



1 



(/ii+2,ii,fei + fi 1 -2,j 1 ,ki 



(47) 



+/ii,ji+2,fei + /ii,ji-2,fci 

+ fn,ji,k 1 +2 + /njiii-2 - 6/ijj^fcJ ■ 

This discrete V 2 is also negative definite and accurate to 
linear order of A. However, with this definition of V 2 , 
if the metric backreaction is negligible and n is an even 
number, a grid point on the lattice will never interact 
with its neighbors. The Fourier-space point of view is 
that k in Eq. (|4"6")l satisfies a stronger PBC: 

jLcff _ jr,eff _ ieff _ ji eff 

^jl+n/2,j 2 ,j3 ~ ft 'jl,i2+n/2j a — rv ji,32,ja+n/2 ~ ^31,32,33 ' 

.(48) 

which indicate that we are only studying modes within 
a smaller fundamental box. The higher modes (interac- 
tion between neighboring points) are not included in this 
discretization scheme. 

Being unsatisfied with the discretization scheme in 
HLATTICE VI. 0, I introduced a new discretization 
scheme in the current version of HLATTICE, HLAT- 
TICE V2.0. The discrete derivatives are defined as 

® l Hi,i2,ia f = 1 o A ^ ifh + iMM ~ 1,»2,*3) 



12A 

-(fil 



+2,i 2 ,i3 



fn-2 



^ 2 k,»2,i 3 f = t 8 (f^^+hia - fh,i2-l,ia) (49) 

~{fil, 12+2, 13 ~ f%l,%2— 2,13)] ) 



8 \ f = - 

3 I»1,»2,*3 ^ _ 12A 



[8 (/ii, 12, 13+1 fii,i2,is— l) 



(/ii, 12, 13+2 fii,i2,is— 2)] 

whose Fourier-space form is 



ik 



iff 

■31,32,33! 31,32,33 1 



fn, 



(50) 



where the effective wave vector kj^ 2 - 3 is 

, eff W • 27Tj! / 2^1 

• 27TJ2 /. 27TJ2 

sin 4 — cos 



sm ■ 



n 

n 



4 — cos 



n 

27TJ3 

n 



(51) 
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TABLE I. Discretization schemes of lattice codes. The "ac- 



curacy" is defined as (V 2 /| 



V 2 / 



Code 



Accuracy 



LATTICEEASY Undefined 


Eq. 


{33} 


0(A a ) 
0(A 2 ) 
0(A 2 ) 


DEFROST Undefined 


Eq. 


cm 


CUDAEasy Undefined 


Eq. 


cm 


HLATTICE V1.0 Eqs. {42} 


Eq. 


(021) 


0(A 2 ) 
0(A 4 ) 


HLATTICE V2.0 Eqs. {49} 


Eq. 


(52} 



This discrete V is two-orders more accurate: V/| d - = 
V/| con , + 0(A 4 ). 

Explicitly written in configuration space, the discrete 
V 2 in HLATTICE V2.0 is 



V 2 / 



/tl+4,*2,*S /il— 4,»2,i3 + /'li 22+4,23 

"T"/ii,i2 — 4,13 ^ fii ,12,13+4 "I" Jii,i2,i3— 4 

"I"l^(/il+l,i2j»3 /il— l,J2,i3 ^~ /il,»2+l>*3 

~^~fh,i2— 1,13 /ii,ia,i3+l + JiiMM—l 
fil+3,t2,i3 /ii— 3,i2,*3 fil,i2+3,i3 
/il, «2— 3,13 /ii,i2,i3+3 f 11,12,13 — ?>) 
+ 64(/ ll+2 ,i 2 ,i 3 + /ii-2,i 2 ,i3 + fii,i 2 +2,i 3 
+ /ii,i 2 — 2,i 3 "I" fii,i2,i3+2 + /ii, 12,13-2) 



-390/ n , 



(52) 



Finally, I summarize all the discretization schemes in 
Table |H 



3. i/oui do uie define the scalar, vector, and tensor on the 
lattice? 



In previous works [74M81L l88l l96l |. the scalar, vector 
and tensor part of /i.y are defined in Fourier space: 



/,;, =-^-[A:f d fcf d -i%(F td ) 2 j A 



where 



i fcf % + i Jfcf d ^ + fcj 



*f d A, = *f d ^ = ^ T = 



(53) 



(54) 



Because A: std does not satisfy the PBC, definition ((53)) 
has to be limited in the fundamental box or in a smaller 
region. It is clear that this definition is not equivalent to 
the discrete form of Eq. {J}, whose Fourier-space counter- 
part is 



where 



+i kfl 3 + i kfA t + hj? , 



k° B Ai — k^hj^ — hj^ — 



(55) 
(56) 



Since the Fourier modes of scalar fields on the lattice 
follow EOM (|3"6"1) . Eq. ([55)1 seems to be a more proper 
definition. Thus, the question is whether one should use 
M TT (k std ) or M TT (k eff ) as the TT projector. One may 
argue that with a reasonable UV cutoff, the difference 
between k and k std is small. However, the TT pro- 
jection is a subtle procedure. Because the scalar metric 
perturbations are typically much larger than the tensor 
ones, we are extracting a small number from a big num- 
ber. A little difference in the projector may lead to a 
big error. With n ~ 10 2 -10 3 and a simple discretization 
scheme, most of the modes in the fundamental box will 
have ~ 1-10% discrepancy between fc std and k . We 
will typically get a numerical GW "noise" with ampli- 
tude A tjno j sc ~ 0.01-0. 1 A s , where A s is the amplitude of 
the scalar metric perturbations. Thus, the relative error 
in the GW amplitude A t is A t , noisc /A t — 0.01/r-0.1/r, 
where r = A t /A s is the tensor-to-scalar ratio of metric 
perturbations. For a model with r < 0.1, we, in princi- 
ple, cannot ignore this potential numerical noise due to 
imperfect TT projection, unless we understand that it 
will vanish in some way. In Sec. IIIIBI this scalar-tensor 
mixing problem due to imperfect TT projection will be 
further discussed with a concrete example. 



C. The sixth-order symplectic-Runge-Kutta hybrid 
integrator 

HLATTICE does not discretize the action along the 
temporal direction. Instead it uses an accurate sixth- 
order symplectic integrator to integrate the EOMs. 

A symplectic integrator is designed to integrate a 
classical system with a Hamiltonian that has the form 
%(p,q) = ^(p) + 'P( c l), where q, p, V(q) and /C(p) are, 
respectively, the general coordinates, the conjugate mo- 
menta, the potential energy and the kinetic energy. I 
have used the curled H to distinguish the Hamiltonian 
from the Hubble parameter H . 

An arbitrary function /(p,q) can be evolved with 



df 

dt 



H/ 



where the functional operator H is defined as 
H/ = {/,«} , 



(57) 



(58) 



with {•,•} being the Poisson bracket. Similarly, we can 
define K = {•, K} and P = {-,V}. Note that H = K + P, 
and that K and P do not commute. 

The solution of (|57| can be formally written as 



f(t + dt) 



„Hdt 



/(*) 



(59) 



The nth-order symplectic integrator is constructed by 



factorizing e 



Hdt _ (K+P)dt 



as 



Hdt _ ciKdt dxPdt c 2 Kdt d 2 Pdt 

o — o o C o ■ 



+ 0(dt n+i ) , (60) 
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where c\,d\, c%,d%, ... are constant c-numbers. 

The operators K and P can be regarded as, respec- 
tively, a Hamiltonian for free particles and that for "in- 
ertialess" particles with a potential. Consequently, the 



exact evolution of the system under 



and 



can 



be achieved numerically. More explicitly, the solutions 



and 



Kdt 



P 



q 



(61) 



(62) 



where §^ 

dp 



should be understood as the vector 

(li"> -) and if is defined in the same way. 

Equations (|6H62p are exact for a finite dt. Using the 
right-hand side of Eq. (|60|) to evolve the system, we will 
only have an error term that scales as dt n+1 . (Strictly 
speaking, there are also machine round-off errors, which 
are ~ 10~ 17 for Fortran double precision numbers that 
are used in HLATTICE.) 

Because symplectic integrators are very stable, they 
are often used to study long-term evolution of many-body 
systems in astronomy and particle physics (97 1 98]. The 
most well-known and oft-used symplectic integrator is 
the second-order one, 



g Hdt _ e Kdt/2 e Pdt e Kdt/2 



0(dt 3 ) , 



(63) 



which is equivalent to the leapfrog algorithm used in 
other lattice codes l&j. lolL l92| . 

In HLATTICE, a modified sixth-order symplectic in- 
tegrator is used. Before introducing the integrator, I 
will write down the discretized Hamiltonian of the scalar 
fields and gravity on the lattice. 

Writing the integrals (I22j|26[) as the sums of the in- 
tegrand on the lattice, we obtain the action of the dis- 
crete system described by (m + 6)n 3 general coordinates, 
<t>l{h,i2,h) and /3ij(ii,i 2 ,i 3 ) 0- < t < m, -n/2 < 
*ii*2,*3 < ^/2), and by their time derivatives. Since 
rescaling the action by a constant factor does not change 
the EOMs of the system, we will drop the factor A 3 in 
the discretized action. 

The conjugate momentum of cj>i(ii, ii\ is) is 



and that of ftj (ii, «2> ^3) is 



(64) 



%J(. 1)i2 , i3 )=^T^ /2 ( 2 -^)(^-^i) 



(i1.i2.i3) 

(65) 



Now we are ready to write down the Hamiltonian of 
the system, given by 



H = /Ci + AC 2 + V 



(66) 



where JC\ is the kinetic energy of the scalar fields and 
the sum of off-diagonal terms in the "kinetic energy" of 
gravity, 



lattice 



/C 2 is the sum of diagonal terms in the "kinetic energy" 
of gravity 



K * = -m £ 



-/3/2 



1 lattice 



2£nL-f£nJ 

i=l \i=l / 



(68) 

and V is the sum of all gradient and potential energy 
terms, 

v = £ eV 2 W{<l>i,<l>2,...,ci> m ) + \g ij d i <t> t d j cl> l 



lattice 



1/3 



\latticc / 

x[£ ($3,1+031,2 + 012,3 
lattice 

— 2/323,1/331,2 — 2^31,2^12,3 — 2^12,3/323,1 

— ^22,1/333,1 — /333, 2 /3ll,2 — /3ll, 3^22,3 

+2^23,2^11,3 + 2031,3022,1 + 2/3i2,1033, 2 )l , (69) 



with g lJ given by Eqs. (|19l) . 

The symplectic integrators found in earlier works |9£ 

Il0l| can not be directly used here. This is due to two 
problems: (i) H = {■, H} contains three noncommutative 
operators Ki = {-,/Ci}, K 2 = {-,/C 2 } and P = {-,V}, 
while in the literature only two-term symplectic factor- 
ization formulas are given; (ii) K 2 is noncanonical, as it 
depends on both /3n and ng„. 

The first problem in principle could be solved by itera- 
tive factorization. We can treat Ki + K 2 as one operator, 
factorize e ( K i+ K 2)dt+Pdt ^ anc | g na v]y factorize each fac- 
tor that contains Ki + K 2 . This procedure, however, is 
not optimal, and leads to a factorization formula with 
hundreds of factors (for sixth or higher order). Indeed 
a much simpler symplectic factorization exists, as I give 
below. 

For arbitrary operators A, B, C, commuting or not, 
e (A+B+c)dt can k e factorized as 



,(A+B+C)<it 



e c 3 Adt/2 e c 3 -Bdt/2 e c 3 Cdt e c 3 &dt/2 e (c 3 +c 2 )Adt/2 

x e c 2 Bdt/2 e c 2 Cdt e c 2 Bdt/2 e (c 2 +c 1 )Adt/2 

x e c 1 Bdt/2 e c 1 Cdt e c 1 Bdt/2 e (c 1 +c )Adt/2 

x e c Bdt/2 e c Cdt e c Bdt/2 e (c +c 1 )Adt/2 

x e ciBdt/2 gCiCrft e ciBdt/2 e ( Cl +c 2 )Adt/2 

x e c 2 Bdt/2 e c 2 Cdt e c 2 Bdt/2 e (c 2 +c 3 )Adt/2 

x £ c 3 Bdt/2 £ c 3 Cdt e c 3 Bdt/2 e c 3 Adt/2 



+0(dt 7 ) , 



(70) 
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where 



ci = -1.17767998417887 
c 2 = 0.235573213359357 , 
c 3 = 0.784513610477560 , 
c = 1 - 2(ci + c 2 + c 3 ) . 



(71) 



Eq. (|70)) can be checked by expanding both sides up 
to sixth-order in dt. A Python script doing this te- 
dious but straightforward job can be downloaded from 
http://www.cita.utoronto.ca/~zqhuang/work/symp6.py . 

The symplectic factorization is not unique. For exam- 
ple, one can interchange A and B in Eq. ([TO)) to generate 
another correct sixth-order symplectic factorization with 
a different 0(dt 7 ) residual term. 

While the proof of a symplectic factorization is always 
trivial, the technique to search for one is not. A general 
approach is to assume a form of the factorization with a 
few undetermined coefficients and solve for these coeffi- 
cients by requiring the exact cancellation of lower-order 
terms. For a high-order symplectic factorization, these 
coefficients often need to be solved numerically, by us- 
ing the Monte-Carlo technique to search for a solution 
in the high-dimensional parameter space. This is how 
I obtained Eq. (JTOJ). If C = 0, Eq. ([70) degrades to 
a two-term symplectic factorization formula containing 
the same set of coefficients in Eq. flTTl) . which, indeed, 
was found two decades ago in Ref. [101| . 

Because factors containing B appear more frequently 
in Eq. ([70)) . I let B = Ki, whose numerical evaluation is 
less expensive. Letting A = K 2 and C = P and noticing 
that the c-numbers can always be absorbed into dt, we 
now only need to write down the explicit algorithm to 
evolve the general coordinates fa and ftij and their con- 
jugate momenta under the operators e Kld *, e K2Ctt , and 
e pdt . For canonical operators e Kldt this is straightfor- 
ward. We take JCi as the Hamiltonian and write down 
the Hamiltonian equations: 



/ 



.Kidt 



(322 

Pas 

/?23 

fa 

Pl2 

Ml Ili,l2,«3 
Hfe \ii,i 2 ,i3 

^-033 lii ,22 ,*3 
^"/^23 l?i ,22 ,i 3 

n A 1 



'«l,l2,i3 

l»i,ia,»3 
lii,z'2,i 3 

I i l ,i2,»3 



31 l*i,*2i»3 



2 Iii,i2,»3 



■ e-^ 3 n^dt)|. . . \ 

V« I lli,l 2 ,»3 1 



* l «l,«2,«3 



/3 23 
An 

/3l2 



322 | 
333 I 



2e~ 



173- ll P- 



rip,, 



2e~ 



n^d* 



£i ( 

2 Ul,l2,l3 

n« 22 + ^dtl . . . 

2 Itl,t2,»3 



n^ 3a + ^dt\ 



23 Ui, 12,13 



P31 Ijj ,i 2 ,13 



Since the quantities used to evolve the system, ft, il^^, 
Tlp 23 , n^gj, II^ 12 and JCi all remain unchanged under 
e Kldt , the above algorithm can be applied for finite dt 
without any ambiguity. For example, one does not need 
to ask whether K,\ on the right-hand side of Eq. (|72[l is 
evaluated at time t or t + dt. 

The same procedure can be applied to derive the 
lattice- version EOMs under the canonical operator e pdt , 
by taking V as the Hamiltonian and writing down the 
Hamiltonian equations. Because of the complexity of G/ 
and G g , the final result is not human-readable. The inter- 
ested reader is referred to the macro files in the released 
HLATTICE package, where the EOMs are defined via a 
series of compiler-preprocessor macros. 

Under the operator e K2dt , the EOMs are 



= K 2 di 



P22 

n /3 2 2 



/Pi! 
P22 

V 



2e~^ 2 m 
M 2 



2e~' 3 / 2 

2^-3/2 
A/ 2 



■(n ft 
n^ 22 

n /3 3 3 



n /3 22 

n &3 



dt 
dt 



(72) 



n fe )dt 

/ 

where the configuration-space label 1^ j 2 ,i 3 is omitted on 
both sides. The rest of the general coordinates and con- 
jugate momenta do not change under e K2dt . 

However, due to the noncanonicality of K 2 , the exact 
solution for Eq. (|73[) with finite dt cannot be achieved. 
This is because JC2 depends on both /?n, /3 22 , ^33 and 
their conjugate momenta. Consequently /3, n^ ll; ILj 22 , 
Tlp 33 and /C 2 are all dynamical on the right-hand side of 
Eq. (f73"[) . Thus an algebraic solution does not exist for a 
finite dt. To solve Eq. ([73"[) I use a fourth-order Runge- 
Kutta integrator with a smaller time step dt' <C dt. This 
Runge-Kutta subintegrator, unlike a global one, does not 
cost extra memory, because the operator K 2 is local (does 
not contain interactions between different grid points). 
Here we are solving n 3 independent sets of ODE, with 
each set containing six coupled ODE. If a global Runge- 
Kutta integrator had been used, the task would then be 
solving 6n 3 all-coupled ODEs, which is numerically much 
more expensive. Since this step is numerically cheap, I 
can make dt' sufficiently small so that the 0(dt' 5 ) error 
from the fourth-order Runge-Kutta integrator does not 
spoil the global 0{dt e ) accuracy. In HLATTICE dt/dt' 
is an adjustable parameter, whose default value is set to 
be 10. 

Both the 0(dt' 5 ) error from the Runge-Kutta subin- 
tegrator and the 0(dt 7 ) term in Eq. ([70[) can be made 
very small by shrinking dt by a factor of a few. This 
symplectic-Runge-Kutta hybrid integrator can thus be 
made very accurate without much additional computa- 
tional cost. In Fig. Q] I show a simulation done on an 
eight-core desktop PC in about half an hour. (All eight 
cores are used, as HLATTICE is an OpenMP parallelized 
code.) The fractional energy noises are suppressed to 
< 10~ 12 . Such noises are ~ 10~ 5 -10~ 3 in other public 
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FIG. 1. Lattice simulation for the preheating model, V = 
+ fsVx 2 , using HLATTICE, where A = KT 13 and 
g 2 /\ = 200. The simulation is initiated at the end of inflation, 
where I define a = 1 and choose the box size L = The 
solid red line is log 10 (_E gra d/i?tot), where i?g r ad is the mean 
gradient energy, and E to t is the mean total energy. The dot- 
dashed black line and dashed cyan line are log 10 (_E pot /i?tot) 
and log 10 (_Ekin/-Etot), where E pot and -Ekin are the mean po- 
tential energy and the mean kinetic energy, respectively. The 



D. The gauge choice problem 

I have written HLATTICE in synchronous gauge for 
practical convenience. In synchronous gauge the gauge 
condition goo = 9m — is local. In other gauges the 
ten metric variables g M „ are constrained by four global 
constraint equations, which have to be solved at every 
time step in order to eliminate the four gauge degrees 
of freedom. This is computationally expensive. Another 
option is to keep all ten variables g^ v on the lattice and 
evolve them with Hamiltonian equations and four exter- 
nal constraint equations defined by the gauge condition. 
However, a symplectic integrator cannot deal with ex- 
ternal constraint equations, at least not in a trivial way. 
Thus, writing HLATTICE in other gauges is a difficult 
problem. Nevertheless, I will discuss the theoretical as- 
pect of a generic gauge choices without implementing it 
in a numerical code. 

At first glance, gauge invariance must be broken when 
we use the exact energy-momentum tensor T^ u on the 
right-hand side of the Einstein equation, while keeping 
only the first-order terms in the Einstein tensor G^ v on 
the left-hand side: 



^i(O) _r_ /^(l) rr-icxact 



(74) 



Here G^} is the background quantity, and contains 



dotted blue line, log 10 \3H 2 M 2 /E tot - 1|, shows that the con- linear terms of metric perturbations Sg^ (or spacetime 



straint equation is satisfied at the 10 level. 



available lattice codes and cannot be suppressed much by 
shrinking dt, because the error term behaves as 0(dt 3 ) in 
those codes. Here, for illustration purposes, the metric 
perturbations are turned off, and a low resolution n — 64 
is chosen. It is about 10 times more expensive to include 
metric perturbations. Hence a simulation with metric 
perturbations typically takes about 5(J|) 3 ( # o/corcs ) ^ 
on a desktop PC, assuming the simulation can be done 
with roughly the same number of (~ 50, 000) time steps. 
However, practically a simulation including metric per- 
turbations is often memory-limited. A simulation with 
n > 256, including metric perturbations, either crashes 
(on low-memory machines) or becomes very slow on most 
computer architectures. This deficiency could be par- 
tially cured by using distributed memory, at the price of 
frequent boundary data exchange. I leave the develop- 
ment of a MPI-parallelized (distributed-memory) version 
of HLATTICE to future work. 



In HLATTICE a few lower-order symplectic integra- 
tors are also given as alternative options. The computa- 
tional cost can be reduced if a lower accuracy is tolerable. 



(21 

derivatives of them). Analogously, we can define G^„, 

(3) 

G/ij/ , etc. In the presence of small-scale nonlinearity in 
T®* act , the spacetime derivatives of the metric perturba- 
tions are not necessarily small, but we still assume 5g^ v 
to be small. (Collapsed objects like blackholes are not 
considered here.) For example, in the late-time Universe 
around a dark matter halo, the Laplacian of the New- 
tonian potential V 2 <J>at can be larger than the "zero-th 
order" quantity H 2 , while <&at remains ~ 10 -5 . Because 
G^J contains d 2 g llv terms, it can be comparable to G^°J. 

(2) 

The second-order G^l , which for dimensional reason does 
not contain more derivatives, is suppressed by one more 
power of 5g^ v . Hence it can be ignored. 

Now, let us perform an infinitesimal coordinate trans- 
formation £ M — > — e£ M . The two sides of Eq. ([74)) 
transform as 



G (0) 

flU 



G« 

^ flV 



G (o) 

flV 



g« 

+ eL*TH 



eL c G^ 



^ftU >( 75 ) 



(76) 



where is the Lie derivative along £. 

If we had followed the (usual) fir st-or der gauge trans- 
formation formulas for the metric [l02| ] . we would have 
discarded the last "second-order" term on the right-hand 
side of (|75[) . However, we know that we should not do 

so, as on small scales eL^G^v is of the same order of 

eL^G ( °}. Both ter ms are needed to match the eL^T^, act 
term in Eq. (|76p. otherwise the Einstein equation in the 
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new gauge would not be equivalent to the one in the old 
gauge. 

In summary, to be self-consistent we need to use new 
gauge transformation rules where spacetime derivatives 
of Sg^ are treated as zeroth-order quantities. This 
then leads to a problem: bg^v themselves may no longer 
be small after a gauge transformation. Physically, this 
means that there are "optimal gauges" where metric per- 
turbations remain small. This is not surprising, as we 
expect, for example, that uniform energy gauge would 
fail in the presence of inhomogeneous matter. 

The question for HLATTICE is then whether the syn- 
chronous gauge is one of the "optimal gauges" . The an- 
swer depends on how long we want to evolve the system 
and how inhomogeneous the scalar fields are. Empirically 
when hij approaches 0(1) we should stop the simulation, 
because in this case we are using rather bad approxi- 
mations. This, however, does not necessarily indicate 
a strong gravitational field. In synchronous gauge one 
has the freedom to choose arbitrary spatial coordinates. 
This corresponds to a nonphysical gauge mode, whose 
amplitude (but not its time derivative) can be eliminated 
at any given moment by a proper gauge transformation 
x i -> x i + C(x). HLATTICE V2.0 includes an option of 
using adaptive spatial coordinates (i.e. eliminating the 
gauge mode in every few steps) . Currently this option is 
labeled as a "beta version" (a trial version) . Its numerical 
stability will be studied in my future work. 

A conjecture is that Newtonian gauge is more optimal, 
since the Newtonian potentials are all physical. I also 
leave the rewriting of HLATTICE in Newtonian gauge 
for future work, if a proper integrator can be found for 
this gauge. 

I end this section by summarizing the current HLAT- 
TICE configuration options in Table fll] 

III. GRAVITY WAVES FROM PREHEATING 

A. Gravity waves from tachyonic preheating after 
hybrid inflation 

In this section I use HLATTICE to calculate gravita- 
tional waves from tachyonic preheating after hybrid infla- 
tion [7l[7i]. Following [Zllzl, I assume a real inflaton 
field <fi and a complex field a = g\ + i 02 . The potential 
reads 

V=\\{a 2 -v 2 ) 2 + l -g 2 ^ 2 a 2 , (77) 

where a 2 = a\ + erf- 

For illustration purpose. I fix the parameters that have 
been used in Figure 4 of [75J: A = 10~ u ,g 2 /\ = 2,v = 

10- 3 V8^M p , and $rU=^ c = 1CT 5 , where <j) c = y/Xv/g is 
the critical point where inflation ends. At the beginning 
of the simulation, the initial metric perturbations are all 
set to zero and a is defined to be 1. The a field is initial- 
ized with random Gaussian fluctuations with "vacuum- 



fluctuation" amplitude Icri^ 2 = |c 2 ,fc| 2 = 1/(2^) and 
l^i.fel 2 = \o- 2 ,k\ 2 = w fc /2, where uj k = ^Jk 2 + m 2 = k, 
as a is massless at the beginning of simulation. This 
is an approximation using the classical lattice simula- 
tion to mimic how these quantum vacuum modes be- 
come classical due to the tachyonic instability. Since the 
growth of da is exponential, this approximation is ex- 
pected to be good. Another problem is that, physically, 
gravity should not respond to these unrenormalized vac- 
uum modes, while on the classical lattice it does respond 
to them. However, since a is light and the lattice UV 
cutoff is not too high, the nonphysical hij excited by 
these initial a fluctuations is negligible. The inflaton <j> 
is set to be initially homogeneous, for two reasons: (i) 
At the beginning, the dominating physics is the tachy- 
onic growth of tr fluctuations. Later <f> fluctuations are 
excited and enhanced by the inhomogenous \ field. The 
vacuum flutuations in cf> remain irrelevant here, (ii) the 
renormalization problem is more severe for the <j) field, as 
at the beginning the energy fluctuations are more sensi- 
tive to <p fluctuations. 

The GW energy spectrum computed using HLATTICE 
is shown in Fig. [2j The fractional energy of GWs per e- 
fold, Sl gw , is defined as 

Sl gw = PaW , (78) 

p cr it d In / 

where / is the GW frequency and p gw is the energy 
density of the GW. Here the critical density, defined as 
p cr it = 3H 2 M 2 , in a spatially flat Universe is the same as 
the mean energy density. I have converted the GW en- 
ergy spectrum to present-day observables in Fig.[2l f2 gWj o 
is defined by replacing all quantities by today's observ- 
ables in Eq. dZSJI. I have used Eqs. (35-36) in Ref. [z| 
to convert the simulation output to present-day observ- 
ables. See also Ref. [79j for a more detailed derivation of 
these formulas. 

In this model, GW are produced during two stages. 

In the first stage, <f> can be approximated as <p — <p c — 
4> c t. The mass square of the a field is approximately 

ml a ~2g 2 cf )c c t . (79) 



The infrared modes k < gy24> c <p c t first start to grow. 
As t increases, more and more modes become tachyonic 
and begin to grow. Bubbles of a field are created in this 
proces s, p roducing gravity waves on roughly the same 
scales (75| . This can be seen from the lower part of Fig. [2] 
In the second stage, <j) condensate is broken due to the 
coupling between <f> and a. The estimation (|79[) is no 
longer valid. The typical scales of inhomogeneity rapidly 
shift towards k ~ g<j>c, producing GW waves on these 
scales. At low k the GW spectrum saturates at a sta- 
tionary level, as shown in the upper part of Fig. [2j The 
saturated low-frequency part of GW is what we are in- 
terested in, as it can potentially be observed with future 
GW probes such as BBO 0,[79|. 
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TABLE II. Available Options in HLATTICE. S2-6 stands for the option of using the second, fourth and sixth symplectic 
integrators; GW stands for the option to calculate gravity waves; r stands for the option of using conformal time (the default 
is physical time); k cff (k std ) stands for the option of using k cff (k std ) to construct the TT projector. 



discretization scheme 

metric 


Eq. (33| 


Eq. 442J) 


Eq. (|49J) 


Minkowski, no perturbations 


S2-6; GW; k std 


S2-6; GW; k ctt ; k std 


S2-6; GW; k ctt ; k std 


FRW, no perturbations 


S2-6; GW; r; k std 


S2-6; GW; r; k ctt ; k std 


S2-6; GW; r; k ctt ; k std 


FRW, with perturbations, synchronous 
gauge, fixed spatial coordinates 


DISABLED 


S2-6; GW; k cft ; k std 


S2-6; GW; k cft ; k std 


FRW, with perturbations, synchronous 
gauge, adaptive spatial coordinates 


DISABLED 


S2-6; GW; k ett 


S2-6; GW; k ett 




5 5.5 6 6.5 

log 10 (f/Hz) 



FIG. 2. Gravity waves from tachyonic preheating after hy- 
brid inflation calculated using HLATTICE VI. with metric 
perturbations and fixed spatial coordinates. The simulation 
resolution is n = 128. The model and parameters are the 
same as those that have been used in Figure 4 of [75}. Here 
figw.o is the fractional energy of GW (per efold in frequency 
/) that would be observed today (assuming radiation domi- 
nation after preheating); h is the current Hubble parameter 
in unit of 100 km s^ 1 Mpc -1 . The outputs are plotted from 
bottom to top per unit time step dt = 0.00313-ff7_j , where 
4> c is the critical point where inflation ends. The lattice sim- 
ulation is initiated at <f> = O.99750 c , where I define a = 1 and 
the box size of simulation L — 0.8H7 ± . 




-30 L lL, i , i , 

1 1.2 1.4 1.6 

a 

FIG. 3. The comparison between energy carried by GW 
and numerical energy noise. The simulation is the same 
one as that described in Fig. O All the quantities are eval- 
uated at the end of the simulation where a m 1.8. The 
solid red line is log 10 (.E grad! ficids/£ ; tot,ficids), where £ gra d,ficids 
is the total gradient energy of </> and a fields, and E to t, fields 
the total energy carried by them. The dot-dashed cyan 
line is log 10 (|G 9 |/Btot, fields), where G g is the "gradient en- 
ergy" of gravity defined in (|26[) . The dashed black line is 
log 10 (_E G w/£ ; to t i fici ds ), where E GW is the energy carried by 
GW. The dotted blue line is \og 10 (H/Et o t, fields), with H be- 
ing the total Hamiltonian given by Eq. (|66[) . 



When gravity is included, the constraint equation is 
always H — regardless of the matter content of the sys- 
tem. This constraint equation can be used to estimate 
the numerical energy noises. In Fig. [3] I plot the en- 
ergy carried by GW and the total Hamiltonian T-L, both 
divided by the energy carried by the scalar fields for com- 
parison. The contribution of energy carried by GW sat- 
urates to ~ 10~ 6 , while the final numerical energy noise 
is r~j 10~ 9 . In other words, energy conservation has been 
checked at the level of about 0.1% of the energy carried 
by GW. 



Comparing the "gradient energy" of gravity G g shown 
in Fig. [3] to the energy carried by GW, we can estimate 
r(k) r~j 1CP 2 for the dominating modes. In this exam- 
ple, the relative contribution of GW due to an imperfect 
TT projector could be as high as ~ 0.1/r, i.e., about 
10 times larger than the physical GW spectrum shown 
in Fig. [5] This provides a possible explanation for why 
the GW spectra found in previous works (see [zl, [93} 
and references therein) are generally larger than what I 
have obtained using HLATTICE. However, the discrep- 
ancy could also be due to the fact that I have included the 
expansion of the Universe, which is ignored in Ref. [75j . 
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Finally, due to the memory limitation, I cannot achieve 
the same high spatial resolution as in Ref. [75| , which the 
authors found is also important for this model. 

The feedback from metric perturbations, however, is 
found to be irrelevant for this model. No significant dif- 
ference has been found in the power spectra of the scalar 
fields between simulations with and without metric per- 
turbations. 



B. Gravity Waves from preheating after chaotic 
inflation 



The example we have shown in the previous section, 
although being observationally interesting, is physically 
complicated. The GW energy spectrum has two peaks 
due to different physics at two stages. Mode-mode cou- 
pling in a wide range of scales becomes important in the 
nonlinear regime, which can hardly be captured by a sim- 
ulation with n = 128. The result shown in Fig. [2] needs 
to be further studied with a simulation with much higher 
resolution if HLATTICE can be MPI-parallelized. 

To better understand the scalar-tensor mixing prob- 
lem, it is better to take a simple example that has been 
well studied and well understood, and can be fully cap- 
tured by a simulation with resolution n — 128. Here I 
take the example of preheating after chaotic inflation, 



(80) 



The parametric-resonance bands and Floquet exponents 
for this model have been studied in detail in Ref. [60(. I 
choose the parameter A = 10 -14 and g 2 /X = 120. This 
set of parameters has been studied in Refs. [Til. [79l. [80| . 
For g 2 /X = 120 the Floquet exponent /z& is shown in 
Figure |4j The dominating mode and the boundary of 
the first resonance band, and the dominating mode in 
the second resonance band are labeled with dashed sky- 
blue, green and cyan lines, respectively. 

I initialize the fields with "vacuum-amplitudes" ran- 
dom Gaussian fluctuations with a cutoff |j2|, |j3| < 
38 (where k off does not significantly differ from k std ). 
This covers the resonance band shown in Figure SJ The 
simulation resolution is n = 128. The box size at the 
beginning of the simulation is L = nA = 20 / (vXfiiw) = 



13.87fl.-,\ where </> ini 



1.714M P is the LATTICEEASY 



default initial background value of <f>, and H ln \ is the ini- 
tial Hubble parameter. The scale factor a is defined to 
be 1 at the beginning of the simulation. Here, in order 
to compare my result with the previous works, I have 
used the same simulation configurations that have been 
used in Refs. 0, III, US] ■ For this model, it is better to 
use conformal time as a time variable. In this simulation 
I thus do not include the backreaction of metric fluctu- 
aions, which requires the coordinate t to be the physical 
time. 

The output of GW energy spectra is shown in Figure [5] 
In the IR part, both the GW projected by M TT (k cff ) 



Hk 0.1 




FIG. 4. The Floquet exponent /ttft(fc) for Xfe( r ) i n the preheat- 
ing model V(4>,x) = j<t> A + ^^X 2 with g 2 /\ = 120, where 
r is the conformal time in units of l/(vA0max) an d k is the 
comoving wave number in units of x/A^max- Here \ = a X an d 
<f> = a<j> are the conformal field values; </f> max is the amplitude 
of background cj> oscillations in the linear regime. The dashed 
sky-blue, green and cyan lines are three characteristic scales 
that will be compared to the simulation results. 



and that by M TT (k std )agree with the results found in 
previous works [zJ, [79|, H(| ■ On the intermediate scales, 
the GW mapped by M TT (k std ) is significantly higher 
than that by M TT (k off ) , though k cff and k std are close 
to each other. In this example, only l%-3% of the 128 3 
wave vectors are in the trustable region, the shaded re- 
gion of the lower panel in Figure [51 where the differ- 
ence between M TT (k std ) and M TT (k cff ) is negligible . 
This is also r dependent. For a problem with smaller 
r, the trustable region can be even smaller. We have 
shown in Sec. IIIBI that at the linear level a mode on 
the lattice exactly follows the EOM with k = k eS and 
that M TT (k cff ) can completely remove the scalar and 
vector components defined by the discrete derivatives. 
Therefore, we may argue that we should trust A/ TT (k cff ) 
rather than M TT (k std ). However, this argument be- 
comes vague in the nonlinear regime. More discussion 
along this line is given in Sec. IIVI Nevertheless, what 
has been explicitly shown here is that caution needs to 
be taken for most of the modes on the lattice. 

In Fig. O the energy spectra of the x ne ld are plot 

"l2I-. |2 , ,.j |2~ 



ted, where n& is defined as = ^ k 2 |XfeT + \x'kY ■ (I 
have used the wavenumber k instead of the frequency ui 
to avoid ambiguity in the nonlinear regime, where w is ill- 
defined.) The growth of in the linear regime (a < 30) 
agrees excellently with the theoretical expectation shown 
in Fig. 21 In the nonlinear regime, the energy spectrum 
still peaks around the resonance band until a ~ 70. En- 
ergy cascading becomes important at a > 70. When the 
energy is peaked in a narrow band, energy cascading is ef- 
ficient. The energy spectrum is soon smoothed at a ~ 80. 
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i ' ' ' ' i ■ ■ ■ ■ i ' ' • 

TT projector using k' 

TT projector using k s 




^^X 2 with A = 10" 14 and g 2 /\ = 120. The simulation 
uses the HLATTICE V2.0 discretization scheme and ignores 
metric feedback. The box size is L = 128A = 20/(V^>ini). 
The blue solid line and red dashed line are GW projected with 

M TT ^efif) and with M TT ^std) ^ respectively In both cages 

a cutoff |ji|, | j7 - 2 1 , |j3 1 < 38 in Fourier space has been used. In 
the upper panel the output at a — 20, 30, 90 are shown 
from the bottom to top. In the lower panel only the outputs 
in the nonlinear regime (from bottom to top, a = 60, 70, 
90) are shown. The shaded regions in the lower panel are the 
regions containing 1% (dark) and 3% (dark and light) of the 
k modes in the fundamental box (—n/2 < ji,j2,js < n/2). 
In both panels the characteristic scales shown in Fig. [4] are 
plotted with the same color code. 



FIG. 6. The energy power spectra of the \ field for the pre- 
heating model V{4>,x) = i^ 4 + with A = 10 ~ 14 and 
g 2 /\ — 120. Here k is the effective wavenumber fc off ; nk is 

defined as n k = ^ \k 2 \ \k | 2 + |Xfe| 2 ] ; (p) is tnc mean energy 
density. The simulation uses the HLATTICE V2.0 discretiza- 
tion scheme and ignores metric feedback. The box size is 
L = 128 A = 2O/(VA0tai). The gray lines from light to dark 
correspond to a = 1, 10, 20, 30, 90. The characteris- 
tic scales shown in Figure [4] are plotted with the same color 
code. 

UV cascading slowly goes on after a ~ 80. However, 
due to the finite resolution of the simulation, UV cascad- 
ing gradually becomes nonphysical on the lattice. When 
the energy spectrum becomes flat, UV cascading will be 
strongly affected by the lattice UV cutoff and should no 
longer be trusted. Physically, UV cascading in a contin- 
uum will continue until it is cut off by the quantum effect 
at very high energy scales where fc 4 is comparable to the 
background energy density. In this sense, the classical 
lattice simulation cannot predict a "final shape" of the 
spectrum. The hope is that, however, other physics such 
as reheating will take place to stop the cascading at some 
diffusion scale. 



IV. DISCUSSION AND CONCLUSIONS 

In Sec. Ill Bl it was shown that a TT projector 
M TT (k std ) leads to a numerical noise in GW, which, 
depending on the tensor-to-scalar ratio and the UV cut- 
off used in the calculation, may or may not be negli- 
gible. This theoretical prediction has been confirmed 
in Sec. MI Bl where the difference between GW ampli- 
tudes calculated with different TT projectors is explicitly 
shown. 

It is also shown that the default TT projector in 
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HLATTICE, M TT (k cff ) , can perfectly remove the scalar 
and vector components defined by the discrete form of 
Eq. (j4j. Since defining the discrete derivatives and 
M TT (k) is a complicated part of the code, it is better 
to check if the actual code agrees with this theoretical 
prediction. In Fig.[7]I show a null test. In the null test, a 
random Gaussian field A with a scale-invariant spectrum 
and r.m.s. amplitude {5k 2 ) 1 / 2 = 1 is generated on a lat- 
tice with resolution n — 128. This field is shown in the 
upper-left panel of each sub-figure. When generating the 
random field A, a cutoff |j 2 |, < 32 is used in sub- 
figure (a2) and (b2), and I.72I, |j3| < 16 in (a3) and 
(b3). The discrete derivatives = didjA are calculated 
using either the HLATTICE VI. discretization scheme 
(left column, sub-figures al-a3) or the HLATTICE V2.0 
discretization scheme (right column, sub- figures bl-b3). 
One component I-F23I is shown in the upper-right panel 
of each sub-figure. The DFT of F^ projected under 
M TT (k std ) are inverse Fourier transformed back to con- 
figuration space. The component IF23I after the TT pro- 
jection is shown in the lower-left panel of each sub-figure. 
Similarly, |F 23 | after TT projection with M TT (k cff ) is 
shown in the lower-right panel of each sub-figure. 

We have confirmed that in the code the scalar compo- 
nent didjA vanishes under the TT projector M TT (k off ). 
(The 1CP 9 level noise is due to the round-off errors in the 
DFT and inverse DFT calculations.) With the same cut- 
off in sub-figures (a3) and (b3), the difference between 
m tt ( k std) and m tt ( k cff) ig much smaller in HLAT- 
TICE V2.0. This again confirms that Vdisc defined by 
Eqs. (|49| is a better approximation to V| con t- 

Although the projector M TT (k off ) perfectly matches 
the configuration-space definition of GWs, caution still 
needs to be taken when we obtain different result us- 
ing M TT (k std ) and M TT (k cff ) . The source of GWs, 
di<t>tdj(j)t, in Fourier space is a convolution, 














FIG. 7. Mapping a test field with TT projectors. Each sub- 
figure contains four panels: the upper-left panel is a 2D slice 
of the test field A; the upper-right panel is a 2D slice of 
I J*23 1 , where Fij = didjA; the lower-left panel is a 2D slice 
of |F 23 | after TT projection under M TT (k std ); the lower- 
right panel is a 2D slice of [.F23I after TT projection under 
M TT (k cff ). In sub-figures (al-a3), the discrete di is defined 
in Eqs. (|42p . In sub- figure (bl-b3), the discrete di is defined 
in Eqs. ((49)) . In sub- figures (a2) and (b2), a Fourier-space cut- 
off U2I, \j3 1 < 32 is applied when generating the A field, 
while this cutoff is | 2 1 , |j3 1 < 16 for sub-figures (a3) and 
(b3). See the text for more details. 



did) ed j 



^3 ^ ' ^ ^ 3xi32^3 



,std 



x (ifcf <i> t ) 



3i '^2 '^3 



stel 



where -n/2 < j[, j' j' 3 , j» , 3% , $ < n/2. The discrete 
Kronecker delta 5^ 3 ^(k) is 1 when k = ^(77,1,772,773) 
(771, 772,773 = 0,±1,±2,...) and zero otherwise. We see 
that while the discrete V is mapped to ik ofr , the mode- 
mode coupling is described by k std in the 5^ function. 
The danger of interpreting k as the physical wave vec- 
tor is that the mode-mode coupling may be inaccurately 
described. 

An ultimate solution to avoid the ambiguity in the TT 
projector might be to evolve everything in Fourier space 
without involving a discrete V approximation. The dif- 
ficulty, however, is that calculating dV/d(j>i in Fourier 
space generally involves convolutions, which is expensive 
for large 77. An attemp t in t his direction is the PSpectRe 
code by Easther et al. [l03j . 

During preheating, the source terms di4>edj(f>£ are com- 
parable to the background energy density. In HLAT- 



2j 3 TfCE the higher-order gravity self-interaction terms < 
Gghij in the Lagrangian are ignored. This is valid at 

(81)least on large scales where the average G g is much smaller 
than the background energy density. On smaller scales 
this might not be a good approximation. In Ref. [93 1 
the authors integrate discretized Einstein equations, and 
find GW enhanced by an order of magnitude when grav- 
ity self-interactions are included. However, the nonlin- 
ear enhancement that they find is on large scales (see 
Figs. 2 and 3 in [HJ]). This is a puzzling result. Note that 
93] suffers from the same scalar-tensor mixing problem. 
Also, as discussed in Sec. IIIB| the discretization of grav- 
ity is not a trivial problem. The authors of [93] find that 
their result is sensitive to the initial conditions of the 
metric. This is a hint that numerical tachyonic instabili- 
ties might exist in their discretization scheme for gravity, 
because physically a weak gravitational field should not 
have chaotic feature. Moreover, numerical noises could 
arise if the integrator is not accurate enough. Ideally 
their results can be checked by adding all the gravity self- 
interaction terms into HLATTICE. However, this will 
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significantly complicate the code and make the simula- 
tions much more expensive. I leave this for future work. 

In HLATTICE the scalar fields are all assumed to be 
canonical. An earlier version of HLATTICE, before met- 
ric perturbations were incorporated, can simulate non- 
canonical scalar fields as well. The noncanonical opera- 
tors are similarly integrated using a Runge-Kutta subin- 
tegrator. I will merge the two versions together in the 
future versions of HLATTICE. The purpose is to accu- 
rately study GW produced in preheating with noncanon- 
ical scalar fields [65| • 

As a general lattice code with a superior accurate in- 
tegrator, HLATTICE can be used in many other fields 
of cosmology. It can be used to study scalar met- 
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